Nonanalytical equation of state of the hard sphere fluid 
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An equation of state of the hard sphere fluid which is not analytical at the freezing density is 
proposed and tested. The nonanalytical term is based on the the classical nucleation theory and 
is able to capture the observed "anomalous increase" of pressure at high densities. It is combined 
with the virial expansion at low densities. 

PACS numbers: 64.10.+h 
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I. INTRODUCTION 

It is well known that liquids can be easily superheated 
or supercooled jj. Equations of state (EOS) describ- 
ing liquids and gases are constructed as analytical func- 
tions for apparent practical and also theoretical reasons 
(thermodynamic approximations usually lead to classical 
mean- field type theories). These functions continuously 
extend to metastable and in most cases also unphysical 
unstable regions. Phase equilibria are then calculated 
from these functions. Yet there are arguments casting 
doubt on this picture. 

The first argument is based on the nucleation pro- 
cess, see Review ^ and references therein. Pure liq- 
uid (any state in general) in a locally stable but globally 
metastable state undergoes after a certain period of time 
a spontaneous transition to a more stable state. In three- 
dimensional systems, liquid/gas and liquid/crystal tran- 
sitions are of the first order and the process is controlled 
by homogeneous nucleation: When a nucleus (droplet) 
of the more stable phase grows above a certain critical 
limit, it keeps growing until bulk liquid freezes or evap- 
orates. Since thermodynamic quantities are time aver- 
ages over functions of configurations, we cannot calcu- 
late them with arbitrary precision. It does not help to 
change the size of the system because in a larger sys- 
tem the probability of the critical fiuctuation is larger, 
the survival time of the metastable state shorter, just to 
compensate higher precision of averages in the larger sys- 
tem. Therefore a value of a thermodynamic quantity in 
a metastable state is not a function at all because it is a 
subject of "essential error" . This error may be viewed as 
a statistical-mechanical analogue of the Heisenberg un- 
certainty principle. 

The second argument is based on the Ising ferromagnet 
which serves as the simplest model of the first-order phase 
transition (in the equivalent form of the lattice g as as a 
model of liquid- vapor equilibrium). It has been proven 
that the magnetization below the Curie temperature 
is not an analytical function of the magnetic field at the 
phase transition point and therefore the magnetization 



cannot be analytically continued to the metastable state. 
(A function is analytical at a point if its Taylor expan- 
sion about this point exists and converges to given func- 
tion in a neighborhood of the point.) Approximate ap- 
proaches to this "essential nonanalyticity" problem are 
based again on the droplet model 0. 

The third hint comes from our experience with the 
hard sphere fiuid. Not only we and other authors Q en- 
countered the impossibility of determining accurately the 
equation of state at large densities (at deeply metastable 
region), but we also found that the data cannot be eas- 
ily fitted to common continuous functions. There is an 
"anomalous" increase of pressure at high densities. In 
contrast, the second derivative of the compressibility fac- 
tor reaches a maximum at high density. 

These arguments led us to an attempt to use formulas 
derived from the simplest version of the classical nucle- 
ation theory in combination with the virial expansion to 
develop an equation of state which is not analytical at 
the freezing point. 



II. THEORY 

A. Classical nucleation theory 

In the classical nucleation theory ^ the Gibbs energy 
of a spherical droplet of solid phase of radius r emerged 
in fiuid is estimated by 



fif) + 47rr^7, 



(1) 



where fj,s and /if are the chemical potentials of the solid 
and fluid phases, respectively, 7 is the s/f interfacial en- 
ergy (surface tension), p is the reduced number density 
(sphere diameter is unity). 

For metastable fluid (p > pfp, where pfp is the fiuid 
density at the freezing point) it holds ps < Pf and the 
Gibbs function G(r) exhibits a maximum at r = r*. 
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where we linearized the difference, {ps — Pf)p ~ —{p — 
Pip)A, A > Q. As soon as the nucleus happens to reach 
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radius r*, it keeps growing. The probability of this spon- 
taneous freezing (per unit time and particle) is propor- 
tional to 



exp 



-AG(r 



kT 



exp 



—const 



(2) 



where const is a positive constant. This is proportional 
to the "metastable uncertainty" , the inherent inaccuracy 
of measurements on the metastable state. 



B. Anomalous behavior in the stable region 

The stable phase (p < ptp) can be viewed as a fluid 
phase with "virtual" nuclei of solid phase of radius r. 
The probability of finding a nucleus of radius r is 



Prob(r) oc exp 



-AG(r) 



(3) 



Quantities as the compressibility factor, Z = p/{pkT) 
(j> denotes pressure, T temperature, and k the Boltz- 
mann constant), are sensitive to the volume of the nu- 
cleus, l^rr^. The anomalous part (caused by "fluctuat- 
ing nuclei") is obtained by integrating over all sizes of 
the nuclei, 



= const' X 
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dr. 



(4) 



Substitution r = a;/(87r7)^/^ and rearrangement of con- 
stants leads to the normalized form 



where a, a > 0, and f3 are constants and 



X exp 



x^\ , 

—t da:. 
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(5) 



(6) 



This integral converges for all t < (then p < pfp in 
Eq. 13, but function ipnit) is not analytical at i = and 
cannot be analytically continued to t > (p > pfp). This 
result is in agreement with similar derivations, e.g., in 0- 
It follows immediately from the above statement that 
the radius of convergence of the virial expansion 
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Bnp" 



(7) 



is less than or at most equal to pfp. In addition, approx- 
imation 

Bn/V""'^ — consto -I- constin + const2n^, 

n > hq, leading to the Carnahan-Starling and other pop- 
ular equations of state Q of the form polynomial(2/) /(I — 
y)'^, is not valid for sufficiently large n and therefore 
such equations cannot describe accurately a vicinity of 
the freezing point. In the above formulas, y = Vp is the 
packing fraction and V is the sphere volume. 
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FIG. 1: Function i/'3(t), Eq. 



C. Beyond stability 

Integral diverges in the metastable region p > pfp 
(t > 0), but it can be approximately evaluated if the up- 
per bound is identified with the maximum r* of G(r). It 
means that configurations with droplets larger than the 
critical size r* are omitted. The error of this approxima- 
tion is estimated by In other words, the nonanalyt- 
icity at p = pfp does not allow extrapolation of Z{p) to 
p > Pfp with precision higher than given by l(2Jl. 

After normalization we get Eq. |SJ with the following 
extension of ||HJ), 
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for t < 
for t > 



X exp 
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Ax. (8) 



This function for rt = 3 is drawn in Fig.^ The maximum, 
appearing also on higher derivatives, matches the third 
hint of the Introduction. 



D. Nonanalytical equation of state 

The equation of state should sew together both the 
low-density region and the anomalous part valid near the 
freezing point. We tried a simple sum of a polynomial in 
p up to order and the anomalous term. 



i=l 



Cp'-^+PMaip-pip)). 



(9) 



where coefficients Ci, i = 1,. ..,to, are determined so 
that the virial coefficients 0,13 up to Em are reproduced. 
Two adjustable parameters a and f3 were fitted to MD 
data j9| up to density pmax- The value of yfp = 7rpfp/6 = 
0.494 was taken from M. The standard deviation a of 
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TABLE I: MD data on the compressibility factor Z at deeply 
metastable states. The error bars are estimated standard de- 
viations in units of the last significant digit, for p > 1.04 
combined with estimated metastable uncertainty. 



p 


Z 


1.02 16.55873(30) 


1.03 


17.2007(20) 


1.04 


17.847(20) 


1.05 


18.59(4) 


1.06 


19.41(8) 


1.07 


20.32(10) 


1.08 


21.30(15) 


1.09 


22.15(30) 



the fit (the objective function) was 

" W 9 7^ ' '^^^^ 

2 — 1 

where ai is the standard error of molecular dynamics 
(MD) datum Zi, see below. 



III. MOLECULAR DYNAMICS DATA 

The values of the compressibility factors at reduced 
densities p < 1.02 are taken from 1^ where also details of 
the molecular dynamics algorithm and finite-size correc- 
tions are explained. The p = 1.03 value was recalculated 
and higher densities newly calculated. We combine runs 
of = 1000, 2000, 4000, and 13500 particles, but at 
lower densities we disregard N = 1000 results because of 
too large finite-size errors (according to Eq. (6) in f£|, the 
influence of periodic boundary conditions is about 0.02 
in Z for the p = 1.04, N = 1000 system). 

The time in MD calculations can be expressed either 
as real time t in reduced time units (mass and sphere 
diameter are unity) or the number of collisions A'coi in the 
system. Without finite size corrections they are related 

by Ea 

7V(3A)i/2(Z~ !)■ 

Two algorithms for obtainingthe initial configuration 
were used. In the first one O"] a periodic cube with 
500 spheres, obtained by Monte Carlo simulation while 
shrinking to the desired density, was replicated. The 
second algorithm was based on soft spheres simulated 
again by Monte Carlo method from random start with 
decreasing temperature until all overlaps disappeared. 
Then during a period of 60-100 time units (about 1000 
collisions per particle) the system was "partially equili- 
brated", i.e., a plateau on the Z{t) convergence profile 
was reached. The second algorithm gave better results 



because in some cases the first one lead to a partly crys- 
taUized cube. 

The metastable MD results need comment because the 
principle inaccuracy limit is reached and data interpre- 
tation is not straightforward. 

At p = 1.02 the system stays at the metastable state 
for 2(1) X 10^" collisions and then it "suddenly" freezes. 
The metastable chunks with approximately constant val- 
ues of Z are long enough to be unambiguously extracted 
and the main source of error is still the statistical uncer- 
tainty of the data; the accuracy could be increased by a 
longer simulation. 

The averaged freezing time for p = 1.03 is 3.3(7) x 
10^ colfisions {t = 1500) for N = 13500 and consis- 
tently 3.8(9) X 10* {t = 5000) for N = 4000 in agree- 
ment with the mechanism of spontaneous homogeneous 
nucleation — the larger system, the higher probability of 
nucleation. However, the N — 1000 system freezes unex- 
pectedly faster, in 6(2) x 10® collisions {t = 1250), which 
we attribute to the influence of the periodic boundary 
conditions. The chunks of data corresponding to the 
metastable fluid are clearly visible, but the points when 
a partly equilibrated metastable state starts and where 
the flrst nucleus of the solid phase appears are not exactly 
deflned. Our analysis of the data using a plot therefore 
contains a subjective factor leading to scattering of 
the Z values by not more than 0.001 if the same proce- 
dure is repeated, exemplifying the principle metastable 
uncertainty. The formally calculated error 0.0014 was 
therefore increased to 0.002. 

For p > 1.04 and partly also for p = 1.04, N > 4000, 
it becomes difficult to determine the metastable parts as 
well as the onset of freezing. There is no plateau on the 
Z{t) curve but rather a shoulder Q. Apparent freez- 
ing appears within time of roughly t = 200 (4000 colli- 
sions per particle). The main source of error in the esti- 
mated Z is the metastable uncertainty which is inevitably 
rather guessed than calculated. Our guesses are more 
pessimistic than those by in addition, our Z data are 
systematically by about 0.001 lower. Our p < 1.04 data 
also match recent results jlQj except the p = 1.04 point 
which is lower, Z = 17.76(2); there are no data in range 
pe [1.05,1.09] in O. 



IV. RESULTS AND DISCUSSION 

In spite of an obscure nature of error estimates in Ta- 
ble we tried to fit them to formula |(2Jl, assuming uni- 
form relative inaccuracies of these errors. The resulting 
estimate of the freezing density, pfp = 0.966(12), is sur- 
prisingly close to the correct value 0.943 8]. 

Results for the nonanalytical EOS © with pmax = 0.95 
(only stable data except the last only slightly metastable 
point) and pmax — 1-09 (all available data) are collected 
in Table ITU It is seen that the fitting is stable and the 
resulting values of a and P are in a physically sensible 
range, especially if we do not consider to < 4 (not enough 
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TABLE II; Parameters a and /3 of the nonanalytical EOS, 
Eq. 0, its residual standard deviation ct, Eq. IIUII . and 
the residual standard deviation (Tan of the expansion in a; = 
y/(l — y) with two adjustable parameters, in dependence on 
the maximum virial coefficient Bm included. 



Pmax =0.95 



Pmax = 1-09 



m 


a 


/3 


a 


0"an 


a 


/3 


a 


fan 


3 


0.6594 


12.81 


347 


36 


0.6341 


13.22 


381 


85 


4 


0.5453 


17.20 


106 


17 


0.5337 


17.69 


106 


33 


5 


0.4740 


23.02 


34 


17 


0.4735 


23.06 


28 


53 


6 


0.4239 


30.76 


11 


4.5 


0.4371 


28.66 


26 


18 


7 


0.3877 


40.56 


5.1 


8.4 


0.4212 


32.31 


31 


25 


8 


0.3602 


52.86 


5.5 


11 


0.4266 


31.06 


30 


22 


9 


0.3407 


66.40 


6.2 


8.6 


0.4619 


23.10 


28 


23 


10 


0.3353 


71.69 


6.3 


21 


0.5473 


11.93 


23 
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N 
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virials to describe low density) and m > 9 ("too stiff" 
virial part), ft is also seen that the equation is not able to 
fit the data within their standard deviations — we would 
need more adjustable parameters. Consequently, using Ui 
in (|10|) may be considered inappropriate. Therefore we 
tried also di = 1 (uniform absolute error) and Ui = Z^, 
(uniform relative error), but the results were qualitatively 
the same. 

A "standard" approach to the HS EOS is a poly- 
nomial va. X = y I (\ — y). The coefficients up to power 
can be determined from the known virial coefficients 
while higher coefficients are fitted to the MD data. Col- 
umn CTan of Table lU contains the residual standard devi- 
ation H1U|) of this analytical function with two adjustable 
parameters. It is seen that both the analytical and non- 
analytical EOSs are for m > 5 comparable. However, 
the result of the nonanalytical EOS are, for m > 7 for 
Pmax = 0.95 or m > 5 for pmax — 1.09, uniformly good (or 
bad) while adding one virial coefficient to the x-expansion 
may worsen the result considerably. 

Similar picture can be obtained from the ability of the 
EOSs fitted in an (almost) stable range, Pmax = 0.95, to 
extrapolate to higher densities. The analytical equation. 
Fig. 121 is better at low densities, but gives nonuniform 
extrapolation to the deeply metastable region, especially 
if many virial coefficients are taken into account. In con- 
trast, the nonanalytical equation, Fig.|31 though worse at 
low densities, gives (for m > 5) uniformly lower values. 
The difference is only a few times larger than the inherent 
metastable error and suggests that there is a systematic 
inaccuracy in the model, in other words, the nonanalyt- 
ical term captures most but not all of the "anomalous 
behavior" . 

Both the analytical and nonanalytical equations ex- 
hibit a "bump" or irregular behavior at p = 1.04. This is 
the density where the well-defined plateaus on the Z(t) 
curves cease to exist and data interpretation contains a 
subjective factor. From the microscopic point of view the 
critical droplet size r* is comparable with the molecule 
size and the nucleation model fails. 



FIG. 2: Analytical equation of state (expansion in a; = y/(l — 
J/)): Deviations of the fit to pmax = 0.95 from MD data. Bm 
is the maximum virial included. Note the change of the scale 
at p = 0.95. 
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FIG. 3: Nonanalytical equation of state 



V. CONCLUDING REMARKS 

In spite of many approximations used (surface energy 
independent on the droplet size, spherical droplets com- 
parable to the atom size, simple combination of the low- 
density expansion and high-density region) the model of 
fluctuating droplets of solid state in fluid is able to pro- 
vide stable results on the equation of state near freezing 
density and enables extrapolation to higher densities. 



Unfortunately, integral i 
for practical application. 



is probably too complicated 
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